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SUMMARY 

This paper presents an extension of a structural finite element mesh 
improvement technique to heat conduction analysis. The mesh improvement con- 
cept was originally presented by Prager in studying tapered, axially loaded 
bars. It was further developed by Kittur, et a 1 . , who showed that an improved 
mesh can be obtained by minimizing the trace of the stiffness matrix. In this 
paper these procedures are extended and applied to the analysis of heat conduc- 
tion in an infinitely long hollow circular cylinder. 


NOMENCLATURE 

N e shape function for element e 

n number of elements 

q temperature gradient, °C/mm 

r radius, mm 

r Q inner radius of the cylinder, mm 

r n outer radius of the cylinder, mm 

r e _i inner radius of element e, mm 

r e inner radius of element e+1 or outer radius of element e, mm 

T 0 temperature at the inner radius, °C 

T n temperature at the outer radius, °C 

j^k? j J stiffness matrix of element e 

k thermal conductivity, W/mm °C 


*Hork funded under NASA Grant NSG-3188. 



thermal conductivity of the material of element e, W/mm °C 
trace of the global stiffness matrix 


Ke 

x 


INTRODUCTION 

The finite element analyst is continually confronted with the decision 
of selecting a "good" mesh for a given problem. Unfortunately there seems to 
be no simple mesh optimization method which is generally applicable. Recently 
Kittur, et a 1 . , (ref. 1) have proposed a method which provides an improved mesh 
over a uniform mesh for a broad class of structural problems. The method is 
based upon the simple procedure of node relocation to minimize the trace of the 
stiffness matrix. Figure 1 presents a flowchart illustrating the procedure. 

An advantage of this method is that the analyst may select the nodal 
positions prior to solving the equilibrium equations. A posteriori methods 
such as those of Shephard (ref. 2) and Carroll (ref. 3) can then be used to 
refine the mesh further. 

Kittur, et al . , have shown that the minimum trace method produces 
excellent results for the problem studied by Prager (ref. 4): a uniformly 

tapered, axially loaded bar. This problem is governed by Poisson's equation 
with Dirichlet boundary conditions. 

In this paper we extend and apply the method to analyze the heat con- 
duction through the cross section of an infinitely long hollow circular cylin- 
der. For this problem, the governing equation is also of the form of Poisson's 
equation. In this case, however, since no heat sources are considered the 
governing equation reduces to Laplace's equation. Both Dirichlet and mixed 
boundary conditions are considered. The results are compared with analytical 
solutions and with finite element solutions obtained with a uniform mesh. 


CONFIGURATION AND PROBLEM DEFINITION 


Consider a hollow circular cylinder of infinite length having inner and 
outer radii: r 0 and r n . Let the thermal conductivity be k, and the temper- 
atures at the inner and outer radii be T 0 and T n respectively. The govern- 
ing equation for the temperature distribution along a radial line is: 


d_ 

dr 
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The boundary conditions are: 

T 


and 
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O 


at 

at 


r 


r 


o’ 


r = r 
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The solution of equation (1) subject to equation (2) is: 


T = 





( 1 ) 


( 2 ) 


(3) 
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Next, suppose that the temperature gradient at the inner surface is 
specified as: q 0 . The boundary conditions are then: 


and 



at r = r 

o 


T = T at r = r 
n n 


( 4 ) 


In this case the solution of equation (1) is: 


T = 


T 


n 


r q 
o H o 



( 5 ) 


FINITE ELEMENT FORMULATION AND MESH OPTIMIZATION 

Figure 2 shows the finite element model. Its consists of a series of 
annular elements. For element (e) let the inner and outer radii be r e and 
r e+ i . The entries of the element stiffness matrix are: 


nr 


[k e .l 

= 2ttk 

Ijj 

e 


r 


e + 1 


e 


dN e dN e 

r -r—!- -r-^ dr 
dr dr 


( 6 ) 


where k is the element conductivity constant and where the element shape 

® p p 

functions and N^ are: 


and 



(7) 


By carrying out the indicated operations the element stiffness matrix 
becomes : 


where S is defined as: 
e 



( 8 ) 


( 9 ) 
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Hence the trace z of the global stiffness matrix is: 


t = 2 


n 


e = l 


where n is the number of elements. 


( 10 ) 


The trace may be minimized with respect to the nodal coordinates by set- 
ting the partial derivative of z with respect to r e equal to zero and solv- 
ing for the ratio r e+ ]/r e : 



( 11 ) 


Since the derivative of a sum is the sum of the derivatives, and assuming that 
the conductivity k is uniform thorugh the body (i.e.: k p _i = k p = k p ,i) 

(11) becomes: 


_a 

9r 


+ r e * r e-l r e+l + r e 

r e " r e-l r e+l ‘ r e 


= 0 


( 12 ) 


which simpl if ies to: 


2r . 2r , 

e-1 e+1 

+ / — = 0 


r e - r e-l / ( r e+l ‘ r e) I 2 
Inverting (13) and factoring: 


(13) 


- 2 e+1 

e 1 r 


r. ? f— - 1 


e- 1 r 


e-1 


e+1 


e-1 


(14) 


Rearranging : 


r e Y r e+1 A = ( r e+\ \ ( r e ^ 


/e-1 A r e 


r e J V r e-1 


(15) 


I 1 I 

Let A = — + and B = — e - , then (15) becomes: 
e e-1 


B( A - l) 2 = A(B - l) 2 


(16) 


or: 


AB( A - B) - (A - B) = 0 


(17) 
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thus : 


(A - B) ( AB - 1) = 0 (18) 

Equation (18) implies that either (A - B) = 0 or (AB - 1) = 0. If (AB - 1) = 0, 

r i r r , 

then AB = 1, which means -^r 1 = 1 = 7 ^ which is impossible. Therefore 

r e r e-l e -1 

(A - B) = 0, which implies that A = B. This leads to the relatively simple 
relation : 


r 


e+ 1 


r 

e 



= Y 


(19) 


where y 

is a constant. The optimal 

nodal 

pos i t ions 

are found through 

repeated 

use of equation (19): 





r 2 r l r 3 r 2 

r e+l 
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e 

r r , 

n n -1 


r l = r 0 = r 2 r l 

r 

e 

r e-l 

r 1 ” r „ 0 ~ Y 

n — 1 n -2 


( 20 ) 


Then: 



( 21 ) 


n n 

— = y Y • - • Y Y - -- Y Y Y=Y 

0 


thus: r n = r QY n &nd Y = Changing n to e, we have. 

r = r nY e or: 
e 0 1 


r = 
e 



( 22 ) 


(23) 


NUMERICAL EXAMPLES 


Example 1 


To illustrate the effectiveness of the method consider an annular cylinder 
with the following temperatures specified on the boundaries and made from a 
material with constant thermal conductivity k = x e = 1.0. 

T = T 0 = 100 °C at r Q = 20 mm 

( 24 ) 


T = T n = 0 °C 


at r n = 50 mm 
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Consider two finite element models, each with four linear elements and 
five nodes: The first model (designated "Uniform Mesh") has equally spaced 

nodes. The second model (designated "Improved Mesh") has its nodal spacing 
governed by equation (23). The objective is the determination of the internal 
temperature distribution. 

The solutions of the finite element governing equations lead to the 
results listed in table I. The error at the nodes is defined as the difference 
between the theoretical results and the finite element results. In order to 
compare the error distribution among the elements, rms errors were calculated 
using 50 uniformly spaced points through the thickness of each element. Each 
element's linear shape function was used to compute the temperature at any 
point along the element from nodal temperature values. Error values were 
calculated at each point from the difference between theoretical and finite 
element temperature values. The rms errors are shown beside the nodal errors 
in table I. The standard deviation of the wi thi n-element rms errors is shown 
at the bottom of table I with the overall rms error for the model which was 
calculated using 200 uniformly-spaced points through the model. (These are not 
quite the same points used for computing the wi thi n-el ement rms errors.) 

The theoretical solution for this example (from eq. 3) is shown in 
figure 3. The nodal values for both models are superimposed on the theoretical 
solution. At this scale, differences between nodal and FEM values are much 
too small to be seen in the figure. The errors for the two models are plotted 
on an expanded scale in figure 4. Since linear elements were used, the errors 
tend to increase towards the center of each element (away from the nodes). 

The uniform mesh has relatively small errors at the nodes and greater 
errors near the inner surface of the annular cylinder where the temperature 
gradient is greatest. The improved mesh is found to have zero errors at the 
nodes and a uniform distribution of errors between elements. The overall rms 
error of the improved mesh is 23 percent less than that of the uniform mesh. 


Example 2 

Next, consider the same cylinder but let the temperature gradient be spec- 
ified on the inner boundary. Specifically, let the boundary conditions be: 


^ = -5.4567833 °C/mm at r = 20 mm 

dr o 

T = T = 0 °C at r = 50 mm 

n n 


(25) 


(The boundary conditions were chosen to yield a nearly identical theoreti- 
cal temperature distribution as in the first example.) 

Table II shows the comparison between the finite element solution and the 
theoretical values of the temperature. Once again the improved mesh shows a 
more uniform distribution of rms errors between elements and the overall rms 
error is lower (by 15 percent) than that of the uniform mesh. The errors for 
the two models are plotted in figure 5. 
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Example 3 


Since the temperature gradient is specified at the inner boundary, it is 
also of interest to know how the values of the temperature gradients obtained 
using the two finite element meshes compare with each other and with the theo- 
retical values, table III provides such a comparison. The errors as plotted 
in figure 6 are much greater (relatively) than in the previous examples, espe- 
cially at the nodes. Note that the largest error is at the inside radius where 
the derivative boundary condition is applied. Generally the Neumann boundary 
conditions are not strongly satisfied in the finite element method (ref. 5). 
Therefore, in general, the FE method is not well suited for gradient-type prob- 
lems. The improved mesh shows a slightly more uniform distribution of rms 
errors between elements and the overall rms error is also slightly lower (by 
7 percent) . 


DISCUSSION 

The numerical examples show that the "improved" mesh defined by equation 
(23) produces results which are in better agreement with the theoretical values 
than those obtained using the uniform mesh. Therefore, the mesh defined by 
equation (12) is an improvement over the uniform mesh for both the temperature- 
fixed boundary conditions and for the mixed boundary conditions. 

In Example 1, the values of the stiffness matrix traces of the uniform and 
improved meshes are 74.6667 and 70.1520, respectively. This indicates that for 
this example the trace is only moderately sensitive to the nodal locations. 

More dramatic differences in the results will occur in problems where the trace 
is more sensitive to changes in the nodal positions. This will occur in prob- 
lems with less regular geometry. 

In the conventional finite element method, if a model with n degrees of 
freedom is refined by adding nodes, a new stiffness matrix needs to be con- 
structed. In general the new stiffness matrix is different in form from the 
original stiffness matrix. However, in the "hierarchical" approach, new ele- 
ments are constructed so that the original n x n stiffness matrix forms a 
submatrix of the new stiffness matrix. In this case, therefore, it is neces- 
sary to compute only the stiffness entries corresponding to the new degrees of 
freedom. This reduces the computing time in assembly. (See Zienkiewicz 
(ref. 6) for a more detailed discussion on hierarchical finite elements.) 

Since error estimation in the hierarchical approach has been established, we 
can take advantage of the error analyses in studying conventional finite 
elements . 


Consider a finite element model with n degrees of freedom (d.o.f). If 
the mesh is to be refined by introducing additional nodes (h-type refinement), 
then it is necessary to know the expected improvement in error before a refine- 
ment step is undertaken. O.C. Zienkiewicz, et al. (ref. 7) and Peano, et al . 
(ref. 8) have shown that if the n+l^h d.o.f. is to be introduced hierarchi- 
cally, then the error in the energy norm is: 


f , - K . u 
, n+1 n + 1 ,n n 


n, 


K 


( 26 ) 


n+1 , n+ 1 
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where f n+ ] is the force corresponding to the n+l th d.o.f., K n+ ] n+1 is the 
stiffness of the n+]th d.o.f <n+i,n is the off-diagonal stiffness relating 
the n+lth d.o.f. to the original n d.o.f. system, and u n is the array of 
nodal displacements of the n d.o.f. system. The subscripts n , 1 of the 
error e refer to the n original d.o.f. and the new d.o.f. 

Zienkiewicz (ref. 6) has used the above error relation to define an error 
indicator in the form: 


n 


2 

n , 1 



^n+1 , n+1 


2 


(27) 


where (, is the finite element residual. 

In an adaptive refinement strategy, these indicators are normally calcula- 
ted for all the d.o.f. corresponding to the next refinement. The indicators 
serve the purpose of identifying the region where refinement is necessary. 

Next, the error corresponding to the previous iteration wherein the n^h 
d.o.f. was added is: 


e 


n-1 


,1 


2 



ill 

K 



The corresponding error indicator is: 


n 


2 

n-1 


,1 





K 


n ,n 


(28) 


(29) 


These derivations are for hierarchical finite elements. However, the 
error with conventional finite elements will be of similar form (ref. 6). 


The most general method of generating a good grid is to have an equal dis- 
tribution of some specified weight function. (See Eiseman (ref. 9) for a com- 
plete discussion on adaptive grid generation.) Often, the error in the finite 
element solution is used as the weight function (ref. 10). Therefore the 
objective is to distribute the error equally among all elements. However, the 
value of the residual C can be obtained only after the equilibrium equations 
are solved. Nevertheless, one way of obtaining an equi-distribution of error 
a priori is by having a uniform element stiffness. As a consequence C will 
be nearly uniform among the elements. The trace minimization procedure deve- 
loped herein produces such a result. Note that each of the ratios in the opti- 
mality condition, equation (19), is a constant y. 


e+1 


n-1 


= Y 


(30) 
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Substituting into equation (9), the element stiffness coefficient is: 


which is a constant. Therefore the trace minimization procedure produces a 
uniform element stiffness. 

Finally, observe the graphs of the errors in figures 4 to 6. As mentioned 
earlier, the errors are more equally distributed with the improved mesh while 
there is a skewed distribution with the uniform mesh. In tables I to III, the 
standard deviation of the wi thi n-element rms errors is lower for the.improved 
mesh than it is for the uniform mesh. Therefore, in all cases, the improved 
mesh distributes the error more uniformly than the uniform mesh. The rms 
errors are exactly equal (within precision of calculation) in the case where 
temperatures are specified at the boundaries. Therefore the mesh obtained in 
this case is optimal. Similar results, however are not obtained in the case 
where both temperature and gradient of temperature are specified because of 
the inability of FEM to strongly satisfy the Neumann boundary conditions. Nev- 
ertheless, it demonstrates the usefulness of the trace minimization procedure 
in a priori grid refinement. 


CONCLUSIONS 

The new method of finite-element grid improvement based on the minimiza- 
tion of the trace of the stiffness method has been extended to the problem of 
heat transfer in a solid body. In elasticity problems, this procedure is 
equivalent to minimizing the potential energy of the model by dividing the 
strain energy equally among the elements. The following conclusions were made: 

1. Nodal positioning obtained by minimizing the trace of the stiffness 
matrix leads to an improved mesh over that obtained by uniform positioning of 
the nodes. 

2. Since trace minimization is an a priori method, the mesh may be refined 
without solving the finite element problem. This makes the minimization proce- 
dure computationally inexpensive to perform. The mesh resulting ftom trace 
minimization may be used as a starting mesh for other mesh refinement procedures 
such as element division or element enhancement (h-methods or p-methods). 

3. The method does not give satisfactory results for analyses in which 
both temperature and gradient of temperature are specified boundary conditions 
due to the inability of the FEM to strongly satisfy Neumann boundary conditions. 
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TABLE I. - COMPARISON BETWEEN UNIFORM MESH, IMPROVED MESH. AND THEORETICAL VALUES FOR 
TEMPERATURE SPECIFIED BOUNDARY CONDITION PROBLEM 


Rad i us , 

Temperature value, °C 

Uni torm mesh error 

Improved 

mesh error 

mm 

Un i form 
mesh 

Improved 

mesh 

Theory 

At node 

In element 
rms 

At node 

In element 
rms 

20. 

25.1487 

27.5 
31.6228 
35. 

39.7635 

42.5 
50. 

100. 

65.3550 

39.0247 

17.7907 

0. 

100. 

75.0000 

50.0000 

25.0000 
0. 

100. 

75. 

65.2453 

50. 

38 . 9260 
25. 

1 7 . 7366 
0. 

0. 

-0. 1096 


0.0000 

0.0000 


1 .04+79 

0.5171 

0.0000 

0.5171 

-0.0988 

-0.0541 

0 . 

0.H>88 

0.4422 

0.2857 

0.0000 

0.5171 

0.0000 

0.5171 

S tandard 

dev i a t i on < 

Df wi th i n-el emen t rms errors 0.-109 

Overall rms error 0.6801 

0.00001 

0.5210 


TABLE II. - COMPARISON BETWEEN UNIFORM MESH, IMPROVED MESH, AND THEORETICAL VALUES FOR 
TEMPERATURE/TEMPERATURE-GRADIENT SPECIFIED BOUNDARY CONDITIONS 


Rad i us , 

Temperature value, °C 

Uni form 

mesh error 

Improved 

mesh error 

mm 

Uni form 
mesh 

Improved 

mesh 

Theory 

At node 

In element 
rms 

At node 

In element 
rms 

20. 

25. 1487 

99.4772 

99.5700 

74.6775 

100. 

75. 

65.2453 

50. 

38.9260 

25. 

17.7 166 
0. 

0.5228 


0.4300 

0.3225 

0.2403 

65.0133 

0.2321 

0.6/86 

27.5 

31.6228 

35. 

39.7635 

49.7850 

24.8925 

0.2150 

0.2980 

38.8207 

0. 1053 

0.4273 

0. 10 76 

0.3782 

17.6977 

0. 

0.0 <89 
0. 

0 . 11)82 
0.2 4 TO 



42.5 

50. 

0. 

0.0000 

0.4694 

Standard 

deviation 

of wi thi n-element rms errors 
Overall rms error 

0 . 1 420 
0 . 4494 

0.0996 

0.3803 


TABLE III. - COMPARISON BETWEEN UNIFORM MESH, IMPROVED MESH, AND THEORETICAL 
TEMPERATURE GRADIENT VALUES FOR TEMPERATURE/TEMPERATURE GRADIENT 
SPECIFIED BOUNDARY CONDITIONS 


Rad i us , 
mm 

Temperature gradient, °C/mm 

Uni form mesh error 

Improved 

mesh error 

Uni form 
mesh 

Improved 

mesh 

Theory 

At node 

In element 
rms 

At node 

In element 
rms 

20. 

25.1487 

27.5 

31 .6228 
35. 

39.7635 

42.5 
50. 

-4.5952 

-3.4923 

-2.8164 

-2.3597 

-2.3597 

-4.8347 

-3.84400 

-3.0578 

-2.4317 

-2.4317 

-5.4568 
-4.3396 
-5.9686 
-3.4512 
-3. 1 182 
-2.7446 
-2.5679 
-2. 1827 

0.8616 

0.4762 

0.3018 

0.2082 

-0.17/0 

0.H029 

0. ’594 

0.2592 
0.11 34 

0.6220 

-0.4947 

-0.3934 

-0.3129 

-0.2490 

0.5183 
0.4122 
0.3278 
0. 1653 

Standard 

dev i a t i on 

of wi th i n-e 1 ement rms errors 0.2081) 

Overall rms error 0.3/55 

0.1490 

0.3478 
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FIGURE 1. - FLOWCHART FOR MESH IMPROVEMENT 
PROCEDURE. 
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ERROR IN TEMPERATURE, °C TEMPERATURE 



RADIAL POSITION, mm 
TEMPERATURE B.C. 


FIGURE 3. - SAMPLE PROBLEM: THEORETICAL AND 
NODAL TEMPERATURE VALUES. 



TEMPERATURE B.C. 

FIGURE 4. - EXAMPLE 1 ( TEMPERATURE/TEMPORARY 
BOUNDARY CONDITIONS): ERRORS IN FEM RESULTS 
FOR CYLINDER TEMPERATURE. 



RADIAL POSITION, MM 
TEMPERATURE/GRADIENT B.C. 



GRADIENT SOLUTION - TEMPERATURE/GRADIENT B.C. 


FIGURE 5. - EXAMPLE 2 (GRADIENT/TEMPERATURE 
BOUNDARY CONDITIONS): ERRORS IN FEM RESULTS 
FOR CYLINDER TEMPERATURE. 


EIGURE 6. - EXAMPLE 2 (6RADIENT/TE MPERATURE 
BOUNDARY CONDITIONS): ERRORS IN FEM RESULTS 
FOR CYLINDER TEMPERATURE GRADIENT . 
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